function [RES,UCoef,SS]=SteadyState_Solver_Res(PP,SS,SS_Res_Setup,Input,UCoef_0)

%% Initiation
if nargin<=3
    if SS_Res_Setup.EqualWageFlag
        Input      =   [1;1;1];
    else
        Input      =   [1;1;1;1];
    end
end

if SS_Res_Setup.EqualWageFlag
    PSI         =   Input(1);
    FracHN_w    =   1;
    DiscFactor  =   Input(2);
    FracHN_wN   =   Input(3);
else
    PSI         =   Input(1);
    FracHN_w    =   Input(2);
    DiscFactor  =   Input(3);
    FracHN_wN   =   Input(4);
end


% PSI         =   1.2;
% FracHN_w    =   1.1;
% DiscFactor  =   0.2;

SS.ExitFlag =   0;
%% Preliminaries
SS.PSI      =   PSI;
SS.BETA     =   1/(1+DiscFactor);
SS.FracHN_w =   FracHN_w;
SS.FracHN_wN=   FracHN_wN;
%--------------------------------------------------------------------------
% Directly Determined by the Steady State Target
%   Note:
%       1. 
%--------------------------------------------------------------------------
% Step 1. Solve mc_p_N, mc_p_H
% SS.mc_p_N   =   ( PP.EPSILON_N-1 )/PP.EPSILON_N + ...
%                 PP.THETA_N/PP.EPSILON_N*SS.Pir_N*(1+SS.Pir_N)*(SS.ir-SS.Pir_N)/(1+SS.ir);
% SS.mc_p_H   =   ( PP.EPSILON_H-1 )/PP.EPSILON_H + ...
%                 PP.THETA_H/PP.EPSILON_H*SS.Pir_H*(1+SS.Pir_H)*(SS.ir-SS.Pir_H)/(1+SS.ir);
SS.mc_p_N   =   ( PP.EPSILON_N-1 )/PP.EPSILON_N;
SS.mc_p_H   =   ( PP.EPSILON_H-1 )/PP.EPSILON_H;
% Step 2. Solve N_H, N_N, w_H, w_N
SS.N_N      =   SS.N*( SS.FracHN_w )/( SS.FracHN_w + SS.FracHN_wN );
SS.N_H      =   SS.N-SS.N_N;
SS.w_N      =   1/( 1+SS.FracHN_wN )/SS.N_N;
SS.w_H      =   SS.FracHN_wN/( 1+SS.FracHN_wN )/SS.N_H;
% Step 3. Solve pY_N, pY_H
SS.pY_N     =   SS.w_N*SS.N_N/(1-PP.ALPHA_N)/SS.mc_p_N;
SS.pY_H     =   SS.w_H*SS.N_H/(1-PP.ALPHA_H)/SS.mc_p_H;

% SS.PAC_N    =   PP.THETA_N/2*SS.Pir_N^2*SS.pY_N;
% SS.PAC_H    =   PP.THETA_H/2*SS.Pir_H^2*SS.pY_H;
SS.PAC_N    =   0;
SS.PAC_H    =   0;
SS.Profit_N =   SS.pY_N-SS.PAC_N-SS.w_N*SS.N_N;
SS.Profit_H =   SS.pY_H-SS.PAC_H-SS.w_H*SS.N_H;
% Step 4. Solve Y, PAC, Profit
SS.Y        =   SS.pY_N+SS.pY_H;
SS.PAC      =   SS.PAC_N+SS.PAC_H;
SS.Profit   =   SS.Profit_N+SS.Profit_H;
% Step 5. Solve B and G
% SS.B        =   SS.Y*PP.B_Y_SS;
% SS.G        =   SS.B/(1+SS.ir)+SS.TAU-SS.B/(1+SS.Pir)-SS.T;

% SS.G        =   0;
% SS.B        =   (SS.G+SS.T-SS.TAU)*(1+SS.ir)*(1+SS.Pir)/(SS.Pir-SS.ir);

% SS.G        =   0;
% SS.B        =   0;
% SS.T        =   SS.TAU-SS.G-SS.B*(SS.ir-SS.Pir)/( (1+SS.Pir)*(1+SS.ir) );

if SS.T+(SS.ir+PP.KAPPA)*PP.b_lb<0
    warning('Transfer is too small to cover the unemployed consumption.');
end

%% Solve the Bellman Equation
% Setup the Approximation Structure

% Initialization
if nargin<=4 || (nargin==5 && isempty(UCoef_0))
    Uhat_0      =   Fun_Vhat(PP,SS,SS.T,SS.Profit,SS.TAU,SS.w_N,SS.w_H,SS.FunApp.U.State);
    UCoef_0     =   SS.FunApp.U.BasMat_0_State\Uhat_0;
end
% Solve Value Function
% BellmanItSetup=   struct('ValItNum',0,'ColItNum',1,'TolItNum',5000,...
%                          'PrintFlag',2,'DampenCoef',1);
ValItNum    =   SS_Res_Setup.Bellman_ValItNum;
ColItNum    =   SS_Res_Setup.Bellman_ColItNum;
TolItNum    =   SS_Res_Setup.Bellman_TolItNum;
DampenCoef  =   SS_Res_Setup.Bellman_DampenCoef;
PrintFlag   =   PP.InfoPrint>=3;
[UCoef,Info,ExitFlag] ...
            =   SteadyState_Bellman(PP,SS,UCoef_0,...
                                    ValItNum,ColItNum,TolItNum,...
                                    PrintFlag,DampenCoef);
if ExitFlag==1
    SS.UCoef        =   UCoef;
    SS.VbarCoef     =   Info.FunCoef.Vbar;
    SS.BellmanInfo  =   Info;
    SS.ExitFlag     =   1;
else
    fprintf('Bellman Solver does not Converge!\n');
    return;
end
%% Solve the Policy
SS.Policy   =   Policy_Solver(PP,SS,SS.T,SS.Profit,SS.TAU,SS.ir_dom,SS.ir_ext,SS.w_N,SS.w_H,...
                              SS.PortfolioInfo,SS.UCoef,SS.PolApp.V.State);
                          
%% Solve the Distribution
% SS.SimuSample=  SteadyState_Distribution_Simulation(PP,SS,5000,100);
% SS.DistApp  =   SteadyState_DistApp(PP,SS);
% Solve the Policy over the Distribution Grid
SS.DistPolicy=  Policy_Solver(PP,SS,SS.T,SS.Profit,SS.TAU,SS.ir_dom,SS.ir_ext,SS.w_N,SS.w_H,...
                              SS.PortfolioInfo,SS.UCoef,SS.DistApp.fzo.State);

% Solve the Approximated Distribution
[SS.TrProb.FlowVec,SS.TrProb.CoefMat,SS.TrProb.UnitTrMat] ...
            =   Distribution_TrProbMat(PP,SS,SS.Pir,SS.dE,SS.PortfolioInfo,SS.VbarCoef,SS.DistPolicy);
SS.Dist_fzo =   SteadyState_Distribution(PP,SS);

%% Aggregation
SS.Agg      =   Aggregation(PP,SS,SS.Dist_fzo.QW_fzo,SS.DistPolicy,SS.TrProb.UnitTrMat);

SS.F_dom    =   SS.Agg.F_dom;
SS.F_ext    =   SS.Agg.F_ext;
SS.Fp_dom   =   SS.Agg.Fp_dom;
SS.Fp_ext   =   SS.Agg.Fp_ext;

SS.B_dom    =   SS.Agg.B_dom;
SS.B_ext    =   SS.Agg.B_ext;

SS.Bp_dom   =   SS.Agg.Bp_dom;
SS.Bp_ext   =   SS.Agg.Bp_ext;

SS.C        =   SS.Agg.C;
SS.C_FI_dom =   SS.Agg.Tot.c.dom;
SS.C_FI_ext =   SS.Agg.Tot.c.ext;
SS.C_RI_N   =   SS.Agg.Tot.c.N;
SS.C_RI_H   =   SS.Agg.Tot.c.H;


SS.L_N      =   SS.Agg.L_N;
SS.L_H      =   SS.Agg.L_H;

SS.FC_Agg   =   SS.Agg.FC;
SS.UI_Agg   =   SS.Agg.UI;

%% Residuals
RES_F       =   (SS.B_dom+SS.B_ext)-PP.B_Total_SS;
RES_L       =   SS.N-(SS.L_N+SS.L_H);
RES_FracHN_L=   (SS.N_H/SS.N_N)-(SS.L_H/SS.L_N);

if SS_Res_Setup.EqualWageFlag
    SS.pC_N     =   SS.pY_N;
    
    SS.FracN_pC =   SS.pC_N/SS.C;
    SS.FracT_pC =   1-SS.FracN_pC;
    SS.FracH_pC =   PP.Cons_FracH*(SS.p_H/SS.p_T)^(1-PP.Cons_ElasHF)*SS.FracT_pC;
    SS.FracF_pC =   (1-PP.Cons_FracH)*(SS.p_F/SS.p_T)^(1-PP.Cons_ElasHF)*SS.FracT_pC;
    
    SS.pC_N     =   SS.FracN_pC*SS.C;
    SS.pC_T     =   SS.FracT_pC*SS.C;
    SS.pC_H     =   SS.FracH_pC*SS.C;
    SS.pC_F     =   SS.FracF_pC*SS.C;

    SS.Cons_FracT=  1-SS.FracN_pC/( (SS.p_N)^(1-PP.Cons_ElasTN) );
    
    RES         =   [RES_F;RES_L;RES_FracHN_L];
else
    SS.FracT_pC =   PP.Cons_FracT*(SS.p_T)^(1-PP.Cons_ElasTN);
    SS.FracN_pC =   (1-PP.Cons_FracT)*(SS.p_N)^(1-PP.Cons_ElasTN);
    SS.FracH_pC =   PP.Cons_FracH*(SS.p_H/SS.p_T)^(1-PP.Cons_ElasHF)*SS.FracT_pC;
    SS.FracF_pC =   (1-PP.Cons_FracH)*(SS.p_F/SS.p_T)^(1-PP.Cons_ElasHF)*SS.FracT_pC;
    
    SS.pC_N     =   SS.FracN_pC*SS.C;
    SS.pC_T     =   SS.FracT_pC*SS.C;
    SS.pC_H     =   SS.FracH_pC*SS.C;
    SS.pC_F     =   SS.FracF_pC*SS.C;
    
    RES_pC_N    =   SS.pC_N-SS.pY_N;
    RES         =   [RES_F;RES_L;RES_FracHN_L;RES_pC_N];
end

RES_FinalGood=  SS.T+SS.Profit+SS.F_ext+SS.F_dom+SS.UI_Agg+(1-SS.TAU)*1 ...
                -SS.C-(SS.B_ext+SS.B_dom+SS.FC_Agg);

SS.RES      =   RES;
SS.RES_FinalGood ...
            =   RES_FinalGood;

%% Solve the Derived Quantities and Prices
% Prices
% SS.p_T      =   (SS.FracTN_pC/(1+SS.FracTN_pC)*1/PP.Cons_FracT)^(1/(1-PP.Cons_ElasTN));
% SS.p_N      =   (1/(1+SS.FracTN_pC)*1/(1-PP.Cons_FracT))^(1/(1-PP.Cons_ElasTN));
% SS.p_H      =   (SS.FracHF_pC/(1+SS.FracTN_pC)*1/PP.Cons_FracH)^(1/(1-PP.Cons_ElasHF))*SS.p_T;
% SS.p_F      =   (1/(1+SS.FracTN_pC)*1/(1-PP.Cons_FracH))^(1/(1-PP.Cons_ElasHF))*SS.p_T;
% Consumption Basket
SS.C_T      =   SS.pC_T/SS.p_T;
SS.C_N      =   SS.pC_N/SS.p_N;
SS.C_H      =   SS.pC_H/SS.p_H;
SS.C_F      =   SS.pC_F/SS.p_F;
% Production
SS.mc_N     =   SS.mc_p_N*SS.p_N;
SS.mc_H     =   SS.mc_p_H*SS.p_H;
SS.Y_N      =   SS.C_N;
SS.Y_H      =   SS.pY_H/SS.p_H;
SS.C_H_s    =   SS.Y_H-SS.C_H;
SS.Z_N      =   SS.Y_N/SS.N_N^(1-PP.ALPHA_N);
SS.Z_H      =   SS.Y_H/SS.N_H^(1-PP.ALPHA_H);
% Government
SS.B        =   SS.Bp_dom;
SS.G        =   SS.B/(1+SS.ir)+SS.TAU-SS.B/(1+SS.Pir)-SS.T;
% Trade Balance
SS.Export   =   (SS.Y_H-SS.C_H)*SS.p_H;
SS.Import   =   SS.C_F*SS.p_F;
SS.NetExport=   SS.Export-SS.Import;
SS.NetExportIdx=SS.NetExport/( SS.Y_H*SS.p_H+SS.Y_N*SS.p_N );

% Capital Flows
SS.NetFinFlow=  SS.B_ext - SS.B_ext*(1+SS.ir_ext)*(1+SS.dE)/(1+SS.Pir);